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Domain decomposition finite element/finite difference 
method for the conductivity reconstruction in a 

hyperbolic equation 

L. Beilina * 


Abstract 

We present domain decomposition finite element/finite difference method for the so¬ 
lution of hyperbolic equation. The domain decomposition is performed such that finite 
elements and hnite differences are used in different subdomains of the computational 
domain: finite difference method is used on the structured part of the computational 
domain and finite elements on the unstructured part of the domain. The main goal of 
this method is to combine flexibility of finite element method and efficiency of a hnite 
difference method. 

An explicit discretization schemes for both methods are constructed such that hnite 
element and hnite difference schemes coincide on the common structured overlapping 
layer between computational subdomains. Then the resulting scheme can be considered 
as a pure hnite element scheme which allows avoid instabilities at the interfaces. 

We illustrate efficiency of the domain decomposition method on the reconstruction 
of the conductivity function in the hyperbolic equation in three dimensions. 


1 Introduction 

With expanding of new compntational technologies and needs of industry it is vital impor¬ 
tance to use efficient computational methods for simulation of partial differential equations 
in two and three-dimensions when computational domains are very large. Domain decompo¬ 
sition methods attracted a lot of interest and is a topic of current research, see, for example, 
[a HD] and references therein. 

It is typical that computational domains in industrial applications often are very large 
and only some part of this domain presents interest. In such cases a domain decomposition 
approach can be attractive when the simple domain discretization can be used in a large 
region and more complex and rehned domain discretization is applied in a smaller part of 
the domain. In this paper we propose to use domain decomposition hnite element/hnite 
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difference approach for the solution of hyperbolic equation which combines flexibility of 
the hnite element method (FEM) and efficiency in terms of speed and memory usage of 
hnite difference method (FDM). To do that we extend a hybrid FEM/FDM method which 
was developed in [5] for the case of acoustic wave equation, to the case of a more general 
hyperbolic equation with two unknown parameters. Similar to our approach in [5], we 
decompose the computational domain such that hnite elements and hnite diherences are used 
in diherent subdomains of this domain: hnite diherence method is used in a simple geometry 
and hnite elements - in the subdomain where we want to get more detailed information about 
structure of this domain. Our goal is to get such a method which will combine hexibility of 
hnite element method and efficiency of a hnite diherence method. It is well known that the 
hnite element method allows to get small features of the structure of the domain through 
the adaptive mesh rehnement. However, this method is quite computationally expensive 
comparing with the hnite diherence method in terms of time and memory usage, see 0-0 
for study of efficiency of these methods. 

In this work we derive explicit schemes for both methods such that hnite element and 
hnite diherence methods coincide on the common structured overlapping layer between these 
domains. Thus, the resulting scheme can be considered as a pure hnite element scheme which 
allows to avoid instabilities at the interfaces [S]. We implement this method in an efficient 
way in the software package WavES [I9] in C++ using PETSc [I6] and MPI (message passing 
interface). 

We illustrate the efficiency of the proposed method on the solution of the hyperbolic co¬ 
efficient inverse problem (CIP) in three dimensions. The goal of our numerical simulations is 
to reconstruct the conductivity function of the hyperbolic equation from single observations 
of the backscattered solution of this equation in space and time. We note, that the domain 
decomposition approach in this case is particularly feasible for implementing of absorbing 
boundary conditions PI- To solve our CIP we minimize the corresponding Tikhonov func¬ 
tional and use Lagrangian approach to do that. This approach is similar to one applied 
recently in [21 El El E] for the solution of different three-dimensional CIPs: we hnd optimality 
conditions which express stationarity of the Lagrangian, involving the solution of a state and 
adjoint equations together with an equation expressing that the gradient of the Lagrangian 
with respect to the conductivity vanishes. Then we construct conjugate gradient algorithm 
and compute the unknown conductivity function in an iterative process by solving in every 
step of this algorithm the state and adjoint hyperbolic equations and updating in this way 
the conductivity function. 

We tested our inverse algorithm on the reconstruction of conductivity function which 
represents small symmetrical inclusions. This problem can be interpreted as the problem of 
the reconstruction of the symmetrical structure of a waveguide and Ending defects in it. Our 
computations show that we can accurately reconstruct large contrast of the conductivity 
function as well as location of all small inclusions using the domain decomposition method 
presented in this work. 

The paper is organized as follows. In Section [2l we present our mathematical model of 
hyperbolic equation and in Section El we describe domain decomposition approach. Energy 
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estimate for the equation of Section|2]is derived in Section 01 In Section0]we formulate state 
and inverse problems for hyperbolic equation. The Tikhonov functional to be minimized and 
the corresponding Lagrangian are presented in Section[ 6 l In Section[7|we describe the domain 
decomposition hnite element/finite difference method to solve the minimization problem of 
Section [ 6 l Finally, in Section [9] we demonstrate efficiency of the domain decomposition 
method on the reconstruction of the conductivity function in three dimensions. 


2 The mathematical model 


The model problem in the domain decomposition method is the following hyperbolic equation 
with the first order absorbing boundary conditions [13] 


1 d^u _ 
dt"^ 

u{x,0) = fo{x), 


V ■ {aVu) = g{x,t), in Qt, 

ut{x,0) = fi{x) in n, 
dnU = —dtu on OVLt) 


( 1 ) 


Here hi C is a convex bounded domain with the boundary dVL G x = (xi, X 2 , X 3 ) G 
and ( 7 ^+“ is Holder space, where fc > 0 and a G (0, 1), c(x),a(x) are the wave speed 
and the conductivity space-dependent functions, respectively. We dehned by fir := x 
(0,T),dflT ■= dfl X (0,T),T > 0. We assume that 

g(x,t) G L2(h2r),/o ^ ^ (2) 


For our purpose we use modihcation of the domain decomposition method developed in 
[S] which was applied for the solution of the coefficient inverse problem for the acoustic wave 
equation in mm- In this work, we use different version of the method used in [H 0] E] , when 
two functions - the wave speed c{x) and the conductivity function a{x) - are introduced in the 
mathematical model of the hyperbolic equation. Similarly with 0-0. we decompose hi into 
two open subregions, flpEM and flpoM such that hi = Qfem^^fdm, and Qfem^^fdm = 0 - 
In flpEM we use finite elements and this domain is such that OVLfem C dflFOM, see hgure 
[2J In VLfdm we will use finite difference method. 

We assume that our coefficients a(x),c(x) of problem ([1]) are such that 

a (x) G [1, di], di = const. > 1, a(x) = 1 for x G Qfdm, 
c (x) G [1, ^ 2 ] , ^2 = const. > 1, c(x) = 1 for x G VLfdm, (3) 

a(x), c(x) G (M^) . 


3 The domain decomposition algorithm 

We now describe the domain decomposition between Qfem and ^fdm domains. This com¬ 
munication is achieved by mesh overlapping across a two-element thick layer around ^Ifem 
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Figure 1: Example of the solution in the domain decomposition between ^fem and ^fdm in one 
dimension with a = c = 1. The interior nodes of the hnite element grid ^fem are denoted by 
stars, while circles and diamonds denote nodes, which are shared between meshes in ^fem and 
^FDM- The circles are interior nodes of the grid in ^fdm, while the diamonds are interior nodes 
of the grid in Q^fem- At each time iteration, solution obtained in ^fdm at uJq is copied to the 
corresponding boundary nodes in ^fem-, while simultaneously the solution obtained in ^fem at 
Wo is copied to the corresponding boundary nodes in ^fdm- 

- see Figure [U First, using the Figured] we observe that the interior nodes of the computa¬ 
tional domain hi belong to either of the following sets: 

ojo Nodes ’o’ - he on the boundary of ^fem and are interior to Qfdm, 

Wo Nodes ’o’ - lie on the boundary of VLfdm and are interior to Qfem, 
w* Nodes ’*’ are interior to ^Ifem, 

UJ+ Nodes ’-f’ are interior to ^fdm 

Then the main loop in time for the explicit schemes which solves the problem ([T]) is as 
follows: 

Algorithm 1 

At every time step k we perform the following operations: 

1. On the structured part of the mesh Qfdm update the FDM solution at nodes w+ and 

Wo- 

2. On the unstructured part of the mesh Qfem update the FEM solution at nodes w* 
and Wq. 

3. Copy FEM solution obtained at nodes Wo as a boundary condition for the FDM solution 
in VLfdm- 

4. Copy FDM solution obtained at nodes Wo as a boundary condition for the FEM solution 
in VLfem- 
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4 Energy estimate 

In this section we prove the uniqueness theorem, or energy estimate, for the function u G 
(fir) of the problem ([T]), using the technique of [T5] . 

Theorem 

Assume that condition ^ on the functions c{x),a{x) hold. Let fl C be a bounded 

domain with the piecewise smooth boundary dQ. For any t G (0,T) let = Q x (0,t). 
Suppose that there exists a solution u G H‘^{VLt) of the problem (Qp. Then the function u 
is unique and there exists a constant A = 74 (||c||q, ||a||Q,f) such that the following energy 
estimate is true for all c, a > 1 in m 



2 

+ \\V^Vu{x,t)\\l^^^^<A 



2 

+ L(o) 


L2(n) 



L2{n) 


( 4 ) 


Proof. 

First we multiply hyperbolic equation in ([T]) by 2dtu and integrate over hi x (0,t) to get 


2 —dttu dtu dxdr — 2V ■ (aVu) dtu dxdr = 2 g dtu dxdr. (5) 




c 

on on 


Integrating in time the first term of ([3]) we get 



0 o 



dt{-^dtU^)dxdT = 


o 2 

on fi o 


diu^ ] {x,t)dx - / -^fpx.t) dx. 


( 6 ) 


Integrating by parts in space the second term of dS]), using conditions Q giving a = 1 
on dVL, and absorbing boundary conditions in ([T]) we get 


V • {aVu) dfudxdr 



0 n 



= 2 / / (dtu) dnudSdr — 2 (aVu) (ydtu) dxdr 



( 7 ) 


0 an 

t 


0 o 



0 an 


(dtu) dSdr — J J adt\VufdxdT. 
0 n 


5 












Integrating last term of ([7]) in time and using initial conditions of the equation ([T]) we 
obtain 



adt\Vu\'^dxdT = — / a\Vu\‘^ {x,t) dx + / a\'Vu\‘^ {x,0) dx 


0 o 


( 8 ) 


= — a\Vu\'^ {x,t) dx + / a\Vfo\‘^ {x) dx. 


Next, collecting estimates ([H]), dZj), ([HD, using the fact that J f (dtu)^ dSdr > 0 and 

0 an 

substituting them in ([HD we have 


dtu^ \ {x,t) dx + / a\Vu\‘^ {x,t) dx 


(9) 



<2 / \g\ \dtu\ dxdr + / —f^{x,t)dx+ / a| V/oP (x) dx. 


0 n 


n 


Finally, to estimate the hrst term in the right hand side of ([HD we use the arithmetic- 
geometric mean inequality 2ab < to obtain 




2 / / Idl ■ \dtu\ dxdr < / / \g\‘^dxdT -|- \dtu\‘^dxdT. 



( 10 ) 


Noting that using ([HD we can write following estimate 

t t 

\dtu\^ dxdr < j j ^—\dtu\‘^ + a\Vdxdr, 
on on 

and substituting the above estimate into fflOli and then the resulting estimate into ([9D we 
have the following estimate 




—dtU^ + a\Vuf]{x,t)dx< / / \g\'^dxdT + i—\dtu\'^ + a\Vu\'^ ] dxdr 



0 n 


0 n 


+ / ^/i+a|V/on 


( 11 ) 


Let us denote 


F{t) := I [ —dtu"^ + a\Vuf ] {x,t) dx. 


( 12 ) 
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We rewrite estimate fllip in the form 


F{t) < A / F(r)(ir + r(t), 


(13) 


where r{t) j J \g\'^dxdT + J + a |V/ot) ix,t)dx. 


0 n 


n 


Applying Gronwall’s inequality to flT3|) with a constant A = A(||c||n, ||a||o,we get 
desired energy estimate ([1]) which also can be written in the form 


n 


-irdtu^ + a I'Vuf ) {x, t)dx < A 


\g\^dxdT+ / {—f^ + a\Vfof){x,t)dx\. 


,0 n 


(14) 


□ 


5 Statement of the forward and inverse problems 

In this section we state the forward and inverse problems. In section [H] we will show how 
these problems can be solved using the domain decomposition algorithm of section [31 

Let the boundary dQ is such that dfl = diQ U 82^ U where Sin and ^2^^ are, 
respectively, front and back sides of the domain hi, and is the union of left, right, top 
and bottom sides of this domain. Let at St ■= difl x (0,T) we will have time-dependent 
observations at the backscattering side of the domain hi. We also define := x 
(0,fi], .Si,2 := d^n X (ti,T), .S2 := <92^ X (0,T) and S3 := dsQ x (0,T). 

We have used 2 model problems in our computations. 

Model Problem 1 

The first model problem is the same as CO but when c{x) = 1 Vx G hi and with non- 
homogeneous initial conditions: 

8‘^ti 

— - V • (aVw) = 0, in Sir, 

u{x,0) = /o(x), ut{x,0) = fi{x) in Q, 

dnU = p(t) on S'!,!, ( 15 ) 

dnU = -dtu on S'!, 2 , 
dnU = -dtU on 5 * 2 , 
dnU = 0 on S3. 

Model Problem 2 
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Our model problem 2 uses homogeneous initial conditions and c{x) = 1 Vx G O in ([T]) 
and is defined as 


— - V ■ (aVu) = 0, in Or, 

M(a:,0) = 0, Mt(x,0) = 0in O, 

dnU = p (t) on Si^i, ( 16 ) 

dnU = -dtu on S'!, 2 , 

dnU = -dtU on 5 * 2 , 

dnU = 0 on S 3 . 

We assume that our coefficient a (x) of problems ffT5|l and ffTHj) is such that 


a (x) E [l,d\, d = const. > 1, a(x) = 1 for x G Qfdm, 
a{x)EC^M.^). 

We consider the following 

Inverse Problem 1 (IPl) Suppose that the coefficient a (x) in the problem satisfies 
Assume that the function a(x) is unknown in the domain VtSflpDM■ Determine the 
function a (x) in m for X G fl\flFDM, assuming that the following function u{x,t) is 
known 

u (x, t) = u {x, t) , V {x, t) G St. (18) 

The question of stability and uniqueness of IPl is addressed in the recent work na. 

Inverse Problem 2 (IP2) Suppose that the coefficient a {x) in the problem /[TSi) satisfies 
Assume that the function a{x) is unknown in the domain Q\Qfdm- Determine the 
function a (x) in m for X G 0\0r_DM, assuming that the following function u{x,t) is 
known 

u (x, t) = u (x, t) , V (x, t) G St. (19) 


6 Optimization method 


In this section we will reformulate our inverse problem IPl as an optimization problem to 
be able to reconstruct the unknown function a(x) in flTHD with best fit to time and space 
domain observations u, measured at a finite number of observation points on Solution 
of IP2 follows from the solution of IPl by taking /o = /i = 0. 

Our goal is to minimize the Tikhonov functional 


F{u, a) 



uY zs{t)dxdt + -7 



OqY dx, 


( 20 ) 


where u is the observed field, u satisfies the equations flT^ . and thus depends on a, Qq is the 
initial guess for a, and 7 is the regularization parameter. Here, zs(t) is a cut-off function. 



which is introduced to ensure that the compatibility conditions at Qt H {t = T} for the 
adjoint problem fl27|) . The function zs can be chosen similarly with [5]. 

For our theoretical investigations we introduce the following spaces of real valued func¬ 
tions 

Hl.= {weH\QT):w{;T)=^0}, 

X X C {Q) , 

f/° = L 2 (fir) X L 2 (nr) X L 2 (fi). 


To solve this minimization problem for model problem 0151) we introduce the Lagrangian 

(aVM)(VA) dxdt 


r/ X X f dXdu ^ ^ 

L{v) = F{u,a) — / ——dxdt + 


' Qj' 


dt dt 


X{x, 0)/i(x) dx 


J VLt 

Xp(t) dadt + 


( 22 ) 


Xdtu dadt -|- / Xdtu dadt, 


' 51,1 


' 5i,2 


'52 


where v = {u,X,a) G U^, and search for a stationary point with respect to v satisfying 
Vh = {u, X, a) G 

L\v,v) = 0, (23) 

where L'{v, ■) is the Jacobian of L at v. 

As usual, we assume that A {x, T) = dtX {x, T) = 0 and impose such conditions on 
the function A that L{u,X,a) := L{v) = F{u,a). We also use the facts that X{x,T) = 
^{x,T) = 0 as well as a = 1 on dO,, together with initial conditions of ffT5|) and boundary 
conditions dnU = 0 on S '3 and 5„A = 0 on fix \ St- The equation fl2^ expresses that for all 
V G U\ 


dL 


dX du 


0 = -;—(«)(A) = — / ——dxdt+ / {a'Vu){VX) dxdt — / X{x,0)fi{x) dx 


dX 




'Si,. 


dt dt 
Xp(t) dadt 




'Si,: 


XdtU dadt + / XdtU dadt, \/XEH\(yLT), 


(24) 


'52 


dL 


dX 


0 = —(m)(m) = / {u — u) u Zs dxdt — / (O)m(x, 0) dx 


du 


' Sx 


’ Qx 


dX du 
dt dt 


dt 


dxdt + / {a'VX){'Vu) dxdt Vm G//^(Ut)- 

J Q.X 


(25) 


Finally, we obtain the equation which expresses that the gradient with respect to a vanish: 

0 = -;—(«) (a) = f CVu)CVX)d dxdt + ^ f (a — ao)ddx, x G fJ. (26) 

da ./o 
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The equation is the weak formulation of the state equation flT^ and the equation 
fl25p is the weak formulation of the following adjoint problem 


d^X 


V ■ (aVA) 


A(-,T) 

dnX 


— {u — u)zs, X e St, 




(•,2^) = o, 


0 , on VLt \ St- 


(27) 


We note that the Lagrangian 0221 ) and the optimality conditions (I2T)) . 02 5 p for the model 
problem 2 will be the same, as for the model problem 1, and only the difference will be that 
these expressions will note have terms with initial conditions. 


7 The domain decomposition FEM/FDM method 

In this section we formulate hnite element and hnite difference methods for the solution of 
model problem 2. FEM for model problem 1 is the same only terms with non-zero initial 
conditions should be induced in the discretization. 


7.1 Finite element discretization 

We discretize SLfeMt = ^fem x (0,T) denoting by Kh = {K} a partition of the domain 
SLfem into tetrahedra K {h = h{x) being a mesh function, defined as h\K = hx, representing 
the local diameter of the elements), and we let Jr be a partition of the time interval (0,T) 
into time intervals J = (tfc-i, ^fc] of uniform length t = — tk-i- We assume also a minimal 

angle condition on the |S]. 

To formulate the hnite element method, we dehne the hnite element spaces Ch, and 
W^. First we introduce the hnite element trial space for u dehned by 

WJi ■-= {we Hi: w\k.j e Pi{K) x Pi(J), VJF G J4, VJ G J.}, 

where Pi{K) and Pi{J) denote the set of linear functions on K and J, respectively. We also 
introduce the hnite element test space dehned by 

:= {we Hi: w\k.j e Pi{K) x Pi(J), VP G P,, VJ G J.}. 

To approximate function a{x) we will use the space of piecewise constant functions Ch C 

L2{n), 

Ch := {u e L,{n) : u\k e Po(P),VP G P,}, (28) 

where Po(P) is the piecewise constant function on P. 

Next, we dehne 14 = x x Ch- The hnite element method now reads: Find Vh G 14, 
such that 

L'{vh){v) = 0Wv eVh- (29) 
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The equation (12 9 p expresses that the hnite element method for the forward problem (1161) 
in flpEM will be: Find Uh G such that VA G and for known ah G Ch, 


dX duh 
dt dt 


dxdt 


dnUhX dxdt + 


lanf 


(a/iVMft)(VA) dxdt = 0, (30) 


and the hnite element method for the adjoint problem (I27P in flpEM reads: Find Xh G W^, 
such that Vh G Wh and for known ah G Ch, 


dX, 

dt 


(O)m(x,0) dx 


'^FEMrj 


dXh du 


dxdt 


(31) 


/ dnXhU dxdt + / {ahVXh){Vu) dxdt = 0. 

'dUpEM J^FEMrp 


We note that usually dim 14 < oo and 14 C as a set and we consider 14 as a discrete 
analogue of the space U^. We introduce the same norm in 14 as the one in 17°, 


ll•||y^ •— ll•ll^ 70 ) (32) 

where Uq is dehned in (l2Tp . From (13^ follows that in hnite dimensional spaces all norms 
are equivalent. This allows us in numerical simulations compute coefhcients in the space Ch- 
However, in the hnite element discretization we write a G 7/2(17) to allow the function a{x) 
be approximated in any other hnite element space. 


7.2 Fully discrete scheme in ^fem 

In this section we present explicit schemes for computations of the solutions of forward and 
adjoint problems in Q^em- After expanding functions Uh{x,t) and Xh{x,t) in terms of the 
standard continuous piecewise linear functions {(pi{x)}f^i in space and {'ipk(t)}h=i in time as 

N M 

EE Uh^,^(Piix)'lpk{t), 

k=l i=l 
N M 

k=l i=l 

where j, and X^ ^ denote unknown coefhcients at the mesh point Xi G Kh and time 
moment t^ G Jr, substitute them into (I5UD and (ITTD . correspondingly, with X{x, t) = u{x, t) = 

We note that we use hnite element method only inside Qeem, and 
thus we will have discrete solutions UhppM '^hppMik ^hpoM ^hpoMik obtained in 
X^FDM as the boundary conditions at dflpEM, after exchange procedure. We obtain the 


Uh{x,t) = 
Xh{x,t) = 
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system of discrete equations: 


N M 


ife+l 


E EE / ^Pi{x)^pj{x) / dti>k{t)dti>i{t) dxdt 

_ 7 . 7 — 1 ^ ^ — 1 K 


K^Q^PEM k,l=l i,j—^ 
N M 


dSdt 


r rtfc+i 

- E EE / dn^i{x)^j{x) / 'lpk{t)'lpl{t) 

dK&dQpEM k,l=li,j=l “'*'=-1 

^ ^ r rtk+i 

+ E EE / ahV 9 ?i(a:)V<^j(a;) / 'ijjk{t)'ipi{t) dxdt = 0. 

For the case of adjoint problem we get the system of discrete equations: 

N M . 


/r^k+1 

iPi{x)ipj{x) / dxdt 

J\^ll.FEM K,l=l 'l,J = L ^ dik -1 

^ ^ P ptk+l 

- E EE ^hi^k / dn^i{x)^j{x) / 'ipk{t)ipi{t) dSdt 

dK<^dnFEMk,l=li ,3 = l “'*'=-1 

^ ^ p ptk+l 

+ E EE ^hi^k ahVipi{x)Vipj{x) j 'ipk{t)i^i{t) dxdt = Q. 

KmFEMk,l=li,j=l 


( 33 ) 


(34) 


Next, we compute explicitly time integrals in (15^ and (154)) using the standard dehnition 
of piecewise-linear functions in time, see [ 1 ] for details of this computation, and get the 
following systems of linear equations: 

M(u^+^ - 2u^ + u^-^) = + -u^ + -u^+^) - + -u^ + -u^+^), 

V ^ I ^6 3 6^ ^6 3 6 ^’ 

M(A^+^ - 2A^ + A^-^) = -r"G(-A^-^ + -A'^ + -A^+^) - + -X^ + -A^+^), 

6 3 6 6 3 6 

(35) 

with initial conditions : 

«(■.«) =f(..0) = 0, (36) 

M-.T) =t(-,r) = 0, (37) 

Here, M is the block mass matrix in space, K is the block stiffness matrix, G is the block 
matrix in space at dVtpEMi is the load vector at time level tk, and A^ denote the nodal 
values of Uh{-,tk) and Xh{-,tk), respectively, r is the time step. 

Now we dehne the mapping Fp for the reference element K such that Fk{K) = K 
and let <,3 be the piecewise linear local basis function on the reference element K such that 
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if o Fx = Then the explicit formulas for the entries in system fl5S]) at each element K can 
be given as: 


^5 = o Fk)k, 

KF = (a/iVv?* o Fk, V^Pj o Fk)k, (38) 

O TV, ‘Pj O FK)dK, 

where {■,-)k denotes the L 2 {K) scalar product and OK is the part of the boundary of element 
K which lies at dVLpEM- 

To obtain an explicit scheme we approximate M with the lumped mass matrix M^, see 
[n] for details about mass lumping procedure. We also approximate terms corresponding 
to the mass matrix in time, and + |A^ + by and A^, 

respectively, which fits to the procedure of mass lumping in time. Next, we multiply fl35|l 
with and get the following explicit method: 

u^+i =(2 - t\M^)-^K)u’‘ - 

^k-i ^ + (2 - 


Finally, for reconstructing a{x) in Qfem we can use a gradient-based method with an 
appropriate initial guess values of oq. The discrete version in space of the gradient with 
respect to coefficient a in fl2B]) take the form: 


9h{x) 


VuhVXhdt + 7(a/i - oo). 


(40) 


Here, Xh and Uh are computed values of the adjoint and forward problems using explicit 
schemes fl39l) . and ah are approximated values of the computed coefficient. 


7.3 Finite difference formulation 


We recall now that from conditions flTT)) it follows that in Qfdm the function a{x) = 1. This 
means that in VLfdm for the model problem flThl) the forward problem will be 


u{x, 0) 


d^u , 
foix), ut{x, 0 ) 

dnU 

dnU 

dnU 

dnU 

dr,.U 


0 , in VLf, 

fi{x) in VLfdm, 
p{t) on Fi,!, 

-dtu on S'i,2, 

-dtu on 5*2, 

0 on S 3 , 

dnUFEM on OSIfEM- 


( 41 ) 
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Then the corresponding adjoint problem in ^fdm will be: 


d^X 


- AA 


A(-,T) 

dnX 

dnX 


—{u — u)zs, X e St, 

0 , on SIt \ St, 
dnXpEM on dVtpEM- 


( 42 ) 


Using standard finite difference discretization of the equation flTT|l in Qfdm we obtain 
the following explicit scheme for the solution of forward problem: 

“5m = - Ufjl, (43) 

and the following explicit scheme for the adjoint problem which we solve backward in time: 


A 


fc-i 

l,j,m 



u)’^ 


l,j,m 


-^AXt^ + 2X'^ 


zs + T AXij „^ 


l,j,m 


\ fc+1 


(44) 


with corresponding boundary conditions for every problem. In equations above, Uij ^ is the 
solution on the time iteration k at the discrete point {l,j,m), {u — is the discrete 

analog of the difference m —ft at the observations points at St, t is the time step, and 
is the discrete Laplacian. In three dimensions, to approximate we get the standard 

seven-point stencil: 


Au 


k 

l,j,m 


“Z+I,j,m - M,j,m + ^ 


2“/ ym + 


u 




6 x 1 

- . +u 


6 Xn 


+ 


l,j,m 


ixl 


(46) 


where 6 xi, 6 x 2 , and 6 x 2 , are the steps of the discrete finite difference meshes in the directions 
xi,X 2 ,X‘i, respectively. 


7.4 Absorbing boundary conditions 

In our domain decomposition method we use first order absorbing boundary conditions [13] 
which are exact for the case of our computational tests of section (9] We note that these 
boundary conditions are implemented in efficient way in the software package WavES [19] in 
the domain decomposition method, and this is the main point of application of this method 
in numerical simulations of section O 

To discretize first order absorbing boundary conditions [13] 

du{x,t) du{x,t) 

-ET = — W- 
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at the outer boundary of ^fdm we use forward finite difference approximation in the middle 
point. This allows to obtain numerical approximation of higher order than ordinary (back¬ 
ward or forward) hnite difference approximation. If we discretize the left boundary of VLfdm 
then we have the condition fl4T]) in the form 


du{x,t) du{x,t) 

dxi dt 

Then the forward hnite difference approximation in the middle point of the above equation 
will be resulted in the following discretization 


u 


fc+i 


U] 


, 3,^11 


+ 


U 


k+1 

7+1, 


U 


/+l,j,m 


U 




U 


l,j,m 


U 


k-\-l 


U 


k-\-l 

l,j,m 


T 


6xi 


Sxi 


= 0, (47) 


which can be transformed to 


^fc+l , k ^ _ k+1 ^ 

“z+1 ,i,m + ^ ^ ^1+1,3,m ^ 

where 5xi is the mesh size in xi direction. For other boundaries of the ^Ifdm absorbing 
boundary conditions can be written similarly. 



7.5 The domain decomposition algorithm to solve forward and 
adjoint problems 

In this section we will present domain decomposition algorithm for the solution of state and 
adjoint equations. We note that because of using explicit domain decomposition FEM/FDM 
method we need to choose time step r such that the whole scheme remains stable. We use 
the stability analysis on the structured meshes im and choose the largest time step in our 
computations accordingly to the CFL stability condition 


T < 



(49) 


Usually, we have the same mesh size 6 xi = 6 x 2 = 6 x 3 = h in {xi,X 2 ,X 3 ) directions, and the 
condition can be rewritten in three dimensions as 


T < h 



(50) 


Algorithm 2 

At every time step k we perform the following operations: 

1. On the structured part of the mesh Qfdm compute from fl45]) . (HID , corre¬ 
spondingly, with and A^, known. 
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2. On the unstructured part of the mesh VLfem compute ^ by using the explicit 

hnite element schemes fl5S]) . correspondingly, with and A^, known. 

3. Use the values of the function at nodes cuo, which are computed using the 

hnite element schemes fl5^ . as a boundary condition for the hnite diherence method 
in VtpDM- 

4. Use the values of the functions X^~^ at nodes Wo, which are computed using the 
hnite diherence schemes (145]) . fl44|) . correspondingly, as a boundary condition for the 
hnite element method in VLpem- 

5. Apply swap of the solutions for the computed functions A^“^ to be able perform 
algorithm on a new time level k. 

8 The algorithm for the solution of an inverse problem 

We use conjugate gradient method for iterative update of approximations a™ of the function 
o/i, where m is the number of iteration in our optimization procedure. We denote 

,r{x) = f VutVXtdt + 7(ar - a„), (51) 

Jo 

where functions Uh {x, t, a™), \h (x, t, a^) are computed by solving the state and adjoint prob¬ 
lems with a := a™. 

Algorithm 3 

Step 0. Choose the mesh Kh in U and time partition of the time interval (0, T). Start with 
the initial approximation a° = Oq and compute the sequences of via the following 
steps: 

Step 1. Compute solutions Uh (x, f, a™) and Xu (x, t, a™) of state ffTHj) and adjoint fl27|) problems 
on Kh and Jr using domain decomposition algorithm of section 17.51 

Step 2. Update the coefficient an := and Jr using the conjugate gradient method 

or' = + aF(x), 

where a, is the step-size in the gradient update HZI and 

(^{x) = -g'^ix) + r(r-^{x), 

with 

^ \\9^-\x)\r 

where dP{x) = —g^{x). 
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o o o o 


a) = flpEM U ^FDM b) flpEM 

Figure 2: a) The hybrid domain = TIfem U ^fdm b) Finite element domain VIfem- 

Step 3. Stop computing a™ and obtain the function ah if either ||5 '™||l 2 (o) < ^ or norms 
IIS'™'! |l 2 (o) are stabilized. Here, 9 is the tolerance in updates m of gradient method. 
Otherwise set m := m + 1 and go to step 1. 

9 Numerical Studies 

In this section we present numerical simulations of the reconstruction of unknown function 
a{x) inside a domain VLfem using the algorithm of section [HI Accordingly to the condition 
f[T7|) this function is known inside ^fdm and is set to be a(a;) = 1. The goal of our numerical 
tests is to reconstruct scatterers of waveguide of figure [2] with c = 4.0 inside every small 
scatterer of figure [21 

The computational geometry hi is split into two geometries, ^fem and VLfdm such that 
n = Qfem U Qfdm, see Figure [2l Next, we introduce dimensionless spatial variables x' = 
x/ (Im) and obtain that the domain Qfem is transformed into dimensionless computational 
domain 

^FEM = {x = (xi, X 2 , X 3 ) e (-3.2, 3.2) X (-0.6, 0.6) x (-0.6, 0.6)} . 

The dimensionless size of our computational domain H for the forward problem is 

H = {x = {xi,X 2 -, X 3 ) G (—3.4, 3.4) X (—0.8, 0.8) x (—0.8, 0.8)} . 

The space mesh in Qfem and in Qfdm consists of tetrahedral and cubes, respectively. We 
choose the mesh size h = 0.1 in our geometries in the hybrid FEM/FDM method, as well as 
in the overlapping regions between FEM and FDM domains. 

In all our computations we use single plane wave p(t) initialized at in time T = [0, 3.0] 
such that 


pit) 


sin (cut), if t G (0, , 

0 , ift>^. 


(52) 


For generation of simulated backscattered data we dehne exact function a{x) = 4 inside 
small scatterers, see Figure [21 and a{x) = 1 at all other points of the computational domain 
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VtpEM- Then we solve the forward problem on a refined mesh which is not the same as used 
in our computations of inverse problem. In a such way we avoid the problem with variational 
crimes. The time step in all our computations is set to be r = 0.006 which satishes the CFL 
condition [ 20 ] . Isosurfaces of the simulated solution for the problem f[T5|l with exact function 
a{x) and cu = 60 in fl5^ are presented in figure O Using this hgure we observe non-zero 
behavior of this solution with initialized initial condition (jSSj)- 

In all our numerical simulations we have considered the additive noise a introduced to 
the simulated boundary data u in flT^ . We have performed following reconstruction tests 
with the same values of parameters in the reconstruction algorithm 3: 

• Test 1. Reconstruction of the function a{x) in Model Problem 1 for uj G [20, 60] in fl52|l 
and additive noise levels a = 3% and a = 10%. 

• Test 2. Reconstruction of the function a{x) in Model Problem 2 for uj G [20, 60] in fl52|) 
and additive noise levels a = 3% and a = 10%. 

In all our tests we start the algorithm 3 with guess values of the parameter a(x) = 1.0 at 
all points in hi. We refer to 0 - 01.00 for a similar choice of initial guess which corresponds 
to starting of the algorithm 3 from the homogeneous domain. The minimal and maximal 
values of the functions a{x) in our computations belongs to the following set of admissible 
parameters 


G {a G U(f2)|l < a(x) < 5}. (53) 

We regularize the solution of the inverse problem by computations with single value of the 
regularization parameter 7 = 0.01 in fl2U]) . Our computational experience have shown that 
such choice of 7 is optimal one in our case. Testing of different techniques, see , for example, 
na, of the computing of regularization parameter is the topic of our ongoing research. The 
tolerance 9 in our algorithm (section [ 8 |) is set to 9 = 10“®. 

We use a post-processing procedure to get hnal images with our reconstructions. This 
procedure is as follows: assume, that functions a^{x) are our reconstructions obtained by 
algorithm 3 of sect ion [ 8 ] where m is the number of iteration when we have stopped to compute 
a{x). Then to get post-processed images, we set 

r a™(x), if a'^(x) > 0.6 max a™(x), 
a™(a:) = < ^fem (54) 

1 1 , otherwise. 

Results of reconstruction for both tests are presented in tables 1,2. Here, computational 
errors in procents are computed for max aj^, where N := m, and are compared with exact 

Ofem 

ones a{x) = 4. 
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a) t = 1.5 


h) t = 1.5 






e) t = 2.1 


f) t = 2.1 
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Figure 3: Test 1. Isosurfaces of the simulated FEM/FDM solution of the problem ()15p with initial 
conditions in ^fem at different times. On a), c), d) we present transmitted data and on b), 
d), f) - backscattered data. 





































































a) Model 1 




e) Models 1 and 2 



d) X 1 X 2 view 
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f) 0:1X2 view 


Figure 4; Test 1. Behavior of noisy backscattered data at time t = 1.8 in both mathematical models 
of section [5j Figure e) presents comparison of backscattered data in both models. 
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Table 1. Results of reconstruction of a(x) for a = 3% together with computational errors in procents. Here, 
N is the final iteration number m in the conjugate gradient method of section\^ 


a = 3% 

a = 3% 

Test 1 

max ajj 

tlpEM 

error, % 

N 

Test 2 

max ajv 

tlpEM 

error, % 

N 

a; = 20 

5.0 

25 

7 

a; = 20 

5.0 

25 

8 

a; = 30 

5 

25 

9 

a; = 30 

4.86 

21.5 

9 

a; = 40 

4.73 

18.3 

10 

a; = 40 

4.19 

4.75 

10 

a; = 50 

5.0 

25 

11 

a; = 50 

5.0 

25 

11 

cu = 60 

5.0 

25 

12 

u! = 60 

5.0 

25 

12 


Table 2. Results of reconstruction of a(x) for a = 10% together with computational errors in procents. Here, 
N is the final iteration number m in the conjugate gradient method of sectionl^ 


a = 10 % 

a = 10 % 

Test 1 

max a^Y 

tipEM 

error, % 

N 

Test 2 

max ajv 

tipEM 

error, % 

N 

a; = 20 

3.29 

17.75 

8 

a; = 20 

3.28 

18 

8 

cu = 30 

4.94 

23.5 

10 

a; = 30 

4.57 

14.25 

10 

cu = 40 

4.35 

8.75 

11 

a; = 40 

3.96 

1 

11 

a; = 50 

4.4 

10 

12 

a; = 50 

4.18 

4.5 

12 

cu = 60 

4.44 

11 

13 

a; = 60 

4.46 

11.5 

13 


9.1 Test 1 

To generate backscattered data at the observation points at Sp in model problem 1, we solve 
the forward problem flT^ . with function p{t) given by 0521) in the time interval t = [0,3.0] 
with the exact values of the parameters a{x) = 4.0 inside scatterers of £gure[2l and a{x) = 1.0 
everywhere else in fl. We initialized initial conditions at backscattered side diQ as 

u{x, 0 ) = exp-(^?+^ 2 +^i) • cost|i=o = exp-(^?+^ 2 +^ 3 ), 
du I 1 2 

— (a:, 0 ) = -exp-(^i+^=+^3) .sint|t=o = 0 . 

Figure [3]-a) presents behavior of this initial condition. 

Figure m presents typical behavior of noisy backscattered data for scatterers of £gure[2]in 
our two models of section [5l Using figure HJ-e) we observe that the difference in the amplitude 
of these two data sets is very small, and thus, we expect that the influence of the non-zero 
initial conditions will not affect to the reconstructions too much. 

Figure |5]-a) presents behavior of the computed norms of differences F'^ = \\5u^ — 
for m > 0, where = \\u^ — u\\zs for a; = 40 in fl52|) and noise level a = 3%. 
Analyzing this hgure for model problem 1 we observe that we achieve convergence in the 
optimization algorithm at iteration m = 10 in the conjugate gradient method. Figures [6] 
presents typical behavior of the computed L 2 norms of differences F™ = — 5u^~^\ \ for 

m > 0 , where = \ \u'^ — u\\zs for different values of uj in fl52]) and noise level a = 10 %. 


21 
























We can see results of reconstruction for model problem 1 in tables 1,2. The reconstructed 
images of the conductivity function for both noise levels and different u are presented in 
figures [71 [9l Figures[7]-b) and|9]-c) show best results of reconstruction which we have obtained 
for cu = 40. We observe that for the noise a = 3% we get correct locations of scatterers and 
values of reconstructed parameter a{x) ~ 4.73 inside them compared with exact one c = 4.0. 
For the noise a = 10% we get correct locations of scatterers and values of reconstructed 
parameter a{x) ~ 4.35 inside them. 

Using these figures we observe that the location of all inclusions in XiX 2 directions is 
imaged very well. However, from hgure [TT] follows that the location in xs direction should 
still be improved. 

9.1.1 Test 2 

This test is similar to the previous one, only we solve IP2 in this case. We start the 
optimization algorithm with guess values of the parameters a{x) = 1.0 at all points in 12 . 
We use the same as in fl5^ set of admissible parameters Ma and the same regularization 
procedure as in test 1 . 

Figure |5]-b) presents behavior of the computed L 2 norms of differences F™ = \\6u‘^ — 
m > 0, for cu = 40 in fl52|) and noise level a = 3%. Analyzing this figure for model 
problem 2 we observe that we achieve convergence in the optimization algorithm at iteration 
m = 10 in the conjugate gradient method. Figure EJ-b) presents typical behavior of the 
computed L 2 norms of differences F™' = \ \5u^ — \ for m > 0 , where (5 m™' = ||m™ — 

for different values of u in (15^ and noise level a = 10 % in this test. 

Results of reconstruction of a(x) for model problem 2 are presented in tables 1,2. Figures 
EM show results of reconstruction for both noise levels and different uj. Using tables 1,2 
we observe that best results of reconstruction are obtained for u = 40. Figures (HJ-b) and 
[Mb) show these results. Using figure | 8 ]-b) we observe that for the noise a = 3% we get 
correct locations of scatterers and values of reconstructed parameter a{x) ~ 4.19 inside 
them compared with exact one a = 4.0 for u = 40. From figure ITOl-b) we see that for the 
noise a = 10 % we get correct locations of scatterers and values of reconstructed parameter 
a{x) ~ 3.96 inside them for u = 40. 

Thus, we again conclude that the location of all inclusions in X 1 X 2 directions is imaged 
very well, but from hgure [TT] follows that the location in X 3 direction should still be improved. 

10 Discussion and Conclusion 

In this work we present domain decomposition FEM/FDM method which is applied for re¬ 
construction of the conductivity function in the hyperbolic equation in three dimensions. 
We have formulated inverse problems and presented Lagrangian approach to solve these 
problems. Explicit schemes for the solution of forward and adjoint problems in the domain 
decomposition approach are also derived. We have formulated different domain decomposi¬ 
tion algorithms: the algorithm 1 describes the overlapping procedure between hnite element 
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a) Model problem 1 


b) Model problem 2 


Figure 5: Differences -F™ = \ \5vJ^ — 6u^ ^|| for w = 40 in (I52p . Noise in backscattered data is 3%. 



a) Model problem 1 


b) Model problem 2 


Figure 6: Behavior of differences = WSvF — dvF ^|| for Model Problem 1 (left figures) and for 
Model Problem 2 (right figures). Noise in backscattered data is 10%. 
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a) cu = 30, max a(x) = 5 

^FEM 


b) cu = 40, max a(x) = 4.73 

^FEM 




c) cj = 50, max a(x) = 5.0 

^FEM 


d) CJ = 60, max a(x) = 5.0 

^FEM 


Figure 7: Test 1. Computed images of reconstructed functions a{x) in model problem 1. We present 
functions d for different u in (1521) and noise level a = 3%. 




a) CJ = 30, max a(x) = 4.86 

^FEM 


b) CJ = 40, max a(x) = 4.19 

^FEM 



Figure 8: Test 2. Computed images of reconstructed functions a{x) in model problem 2. We present 
functions o for different cj in (j52l) and noise level a = 3%. Here we have initialized zero boundary 
conditions in the generation of backscattered data and in the optimization algorithm. 
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a) u! = 20, max a(x) = 3.29 

Qfem 


b) cj = 30, max a(x) = 4.94 

^FEM 




c) o; = 40, max a(x) = 4.35 

flFEM 


d) o; = 50, max a(x) = 4.4 

^FEM 


Figure 9: Test 1. Computed images of reconstructed functions a{x) in model problem 1. We present 
functions d for different u in (j52h and noise level a = 10%. 


P P ^ ^ pi 


a) a; = 30, max a(x) = 4.57 

^FEM 



b) o; = 40, max a{x) = 3.96 

^FEM 




c) u! = 50, max a(x) = 4.18 

^FEM 


d) o; = 60, max a{x) = 4.46 

^FEM 


Figure 10: Test 2. Computed images of reconstructed functions a{x) in model problem 2. We 
present functions d for different uj in (1521) and noise level a = 10%. 
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e) Test 1: cu = 40, a = 10% 



f) Test 2: u = 40, a = 3% 


Figure 11: Computed images of reconstructed functions a(x) in both model problems. We present 
functions a for u; = 40 in ()52p . We observe that reconstruction in x^, direction should be improved. 


and finite difference domains, the algorithm 2 presents solntion of the forward and adjoint 
problems using the domain decomposition FEM/FDM methods, and the algorithm 3 de¬ 
scribes the conjugate gradient algorithm for reconstruction of the conductivity function with 
usage of algorithms 1,2. 

In our numerical tests we have obtained stable and good reconstruction of the conductiv¬ 
ity function a{x) in the range of frequencies oj G [20, 60] . Using tables 1,2 we can conclude 
that the best reconstruction results are obtained in model problem 2 for cj = 40. We can 
also conclude that in all tests we could reconstruct size on xiX 2 -directions for a, however, 
size in x^ direction should be still improved. Similarly with mum we plan to apply an 
adaptive hnite element method in order to get better shapes and sizes of the inclusions in 
all directions. 
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